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Abstract. We present results of global 3D MHD simulations of disk accretion 
to a rotating star with dipole and more complex magnetic fields using a Godunov- 
typc code based on the "cubed sphere" grid developed earlier in our group. We 
describe the code- and the grid and show examples of simulation results. 



1. Introduction 

A wide range of stars have significant intrinsic magnetic fields. Stars are usually 
strongly magnetized during the protostellar stage (T Tauri stars), and after 
collapse to a white dwarf or a neutron star. Many observational properties of 
these stars are determined by the interaction of the accreting disk matter with 
the rotating magnetosphere of the star (see e.g., Bouvier et al. 2007 for review). 
In general, the magnetic axis of the star does not coincide with the rotational 
axis, due to which the magnetospheric flow is complicated and the problem 
requires global 3D MHD simulations. In addition, the magnetic field of the star 
may have a complex structure, which adds complications and the necessity to 
consider this problem in a global MHD approach. 

To solve this problem we developed a special 3D MHD code on the cubed 
sphere grid (Koldoba et al. 2002, see also Putman & Lin 2007) which is some- 
what similar to Yin- Yang grid (Kageyama and Sato 2004). This grid has a 
number of advantages over spherical or Cartesian grids. A "cubed sphere" grid 
had been originally developed for the surface of a sphere for geophysical appli- 
cations (Sadourny 1972; Ronchi, Iacono, & Paolucci 1996). In contrast with 
these authors, we perform simulations in three-dimensional space. We used a 
Godunov-type numerical scheme (Powell et al. 1999; Kulikowskii, Pogorelov, & 
Semenov 2001) and were able to perform pioneering simulations of disk accretion 
to magnetized stars with inclined dipole geometry (Romanova et al. 2003, 2004; 
Kulkarni & Romanova 2005). In this paper we show more recent simulation 
results obtained with our "cubed sphere" grid. 
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Figure 1. Top: cubed sphere grid used in simulations. The grid consists of 
6 blocks corresponding to the 6 sides of a cube. Bottom: Cubed sphere grid 
with high resolution near the equator. 

2. Numerical Method and "Cubed Sphere" Grid 

We consider disk accretion to a rotating magnetized star. This problem is diffi- 
cult to treat numerically because the magnetic field varies strongly with distance 
from the star (~ 1/B? in case of the dipole field and even more steeply for higher 
multipoles), and it is rapidly varying in the laboratory inertial reference frame. 
To minimize errors in calculating the magnetic force, the magnetic field B is de- 
composed into the "main" dipole component of the star, Bo, and the component 
Bi induced by currents in the disk and in the corona (Tanaka 1994). Another 
difficulty with this problem is that the dipole moment changes with time. It 
rotates with angular velocity fi so that the "main" field Bo also changes with 
time. Consequently, in the induction equation there is a large term involving 
Bo . To overcome this difficulty we use a coordinate system rotating with angular 
velocity fi: 

dp/dt + V • (jjv) = 0, 
d(pv) /dt + V • T = pg + 2p v x n - p n x (ft x R), 
d{ P S)/dt + V • (pSv) = , 
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dB/dt = Vx(v x B), 

where v is velocity of plasma in the rotating frame, B is the magnetic field, 
and S is the specific entropy. T is the stress tensor with components = 
P$ik + pviVk + (B 2 5ik/2 — BiB^/Air + Tjfc. Here is viscous stress, p is the 
gas pressure. We consider that the viscous stress is determined mainly by the 
gradient of the angular velocity because the azimuthal velocity is the dominant 
component in the disk. We use the a-viscosity model of Shakura and Sunyaev 
(1973) with the coefficient of dynamic viscosity rjt = ctp/Qx, where a is a 
dimensionless coefficient, a < 1. The viscosity acts only in the accretion disk so 
that the dominant contribution to the viscous stress arises from the gradient of 
the azimuthal velocity (approximately Keplerian) of the plasma. In cylindrical 
coordinates the non-zero components of the viscous stress tensor are 

dn K 3 

T r <j> = = -Vi-faT = IjP , 

where VLk is the Keplerian angular velocity at the given location. We calculate 
momentum fluxes due to the viscous stress at faces of the grid after transforming 
to the Cartesian coordinates. 

The "Cubed Sphere" Grid. The three dimensional grid consists of a set of con- 
centric spheres of radii Rj in a geometric progression with j = 1..Nr. The grid 
on the surface of the sphere consists of six sectors with the grid on each sector 
topologically equivalent to the equidistant grid on the face of a cube. In each 
sector the grid of N x N cells is formed by the arcs of great circles separated by 
equal angles. This grid gives high spatial resolution close to the star which is 
important for our study. Recently we incorporated the option to allow the grid 
to be compressed towards the equatorial plane, to have higher grid resolution 
in the disk (Fig. 1, bottom plots). Such a grid is needed when higher reso- 
lution is required in the disk. Typical grid resolutions used in our simulations 
vary from Nr x A^ 2 =72 x 31 2 cells in each of the six sectors, up to Nr x iV 2 
=288 x 121 2 depending on the problem. The cubed sphere grid naturally lends 
itself to division into 6 x N regions (with N cuts in the radial direction), which 
are calculated in parallel using from 48 up to 240 processors. 

Godunov-Type Finite-Difference Scheme. All variables are evaluated at the cen- 
ters of the cells, and all vector variables are expressed in terms of their Cartesian 
components. Finite difference equations are written for the Cartesian compo- 
nents of vector variables. The finite difference scheme of Godunov's type has 
the form: 

^ V + s m,Tm = Q ■ 

m=1..6 

Here, U = {p, pv,B, pS} is the "vector" of the densities of conserved variables; 
T m is the "vector" of flux densities normal to the face "m" of the grid cell, 
s m is the area of the face "m", V is the volume of the cell, Q is the intensity 
of sources in the cell, and At is the time step. To calculate the flux densities 
J- m , an approximate Riemann solver is used, analogous to the one described by 
Powell et al. (1999) and by Kulikovskii et al. (2001). 
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3. Accretion in stable and unstable regimes 

One of the most striking recent results obtained in our group is the discovery 
of accretion through the interchange instability (Kulkarni & Romanova 2008; 
Romanova et al. 2008). Simulations have shown that a magnetized star may be 
either in the stable or unstable regime of accretion. In the stable regime matter 
accretes to the star in two ordered funnel streams and produces ordered hot 
spots on the surface of the star, leading to periodic light curves (see Fig. 2, left 
panels). In the unstable regime matter penetrates through the magnetosphere 
due to the interchange instability forming a small number (2 to 7) of "tongues" 
(see also Li &; Narayan 2004) which form chaotically at different parts of the 
inner disk, and the light-curve from the resulting stochastic hot spots is expected 
to be irregular (see Fig. 2, right panels). 

There are a number of factors which determine the regime of accretion. If 
the magnetic axis is inclined with the rotation axis at a large enough angle, 
> 30°, then the flow is usually stable (Kulkarni &; Romanova 2009). If the 
inclination is small, e.g., = 5° (used in the majority of our simulations), then 
the stability of accretion depends on various other factors. 

We compared our simulations with a few relevant theoretical approaches. 
The basic theory states that a homogeneous vertical field at the disk-magnetosphere 
boundary does not damp azimuthal perturbations, and therefore is not an obsta- 
cle for the development of unstable </>-modes (e.g. Arons & Lea 1976). A more 
general criterion for magnetized accretion disks states that the disk is unstable 
to growth of </>-modes if 7^ s = —g e ffd\n(Yi/B z )/dr > 2(r m dd/dr) 2 = 7^ (e.g. 
Spruit, Stehle & Papaloizou 1995; see also Kaisig, Tajima & Lovelace 1992). 
Here, —g e ff = r {^K ~ * s the effective gravity, and S = 2ph is the surface 
density. That is, for the instability to start, the surface density per unit mag- 
netic field strength T,/B z should drop off fast enough in the direction of the star, 
that the term 7^ is larger than the term associated with the shear, 7^, which 
tends to suppress the instability by smearing out the perturbations. In addi- 
tion the term —g e ff should be large enough and positive for the instability to 
start. In one set of runs we fixed rotation rate of the star and the initial surface 
density in the disk and varied the accretion rate by varying the a-parameter: 
M ~ £17 ~ a. We observed that cases with a < 0.04 (small M) correspond 
to stable accretion, while cases with a > 0.04 (large M) correspond to unstable 
accretion. We observed that at larger M, the gradient dhx{Yi/B z )/dr is larger 
and the shear 2(r m d£l/dr) 2 is smaller. In addition, at larger M the inner disk 
comes closer to the star, which increases —g e ff- All these factors make instabil- 
ity more favorable. Thus, we observed that increasing in accretion rate leads to 
transition from the stable to the unstable regime (Kulkarni & Romanova 2008; 
Romanova et al. 2008). In another set of runs we fixed the accretion rate in the 
disk by fixing a at some small value, a = 0.02, but varied the rotation rate of 
the star f2. We observed that at low stellar rotation rates the flow is unstable, 
while at higher rates it becomes stable. This is ascribed to the decrease in the 
effective gravity, —g c s, by fast rotation. The accretion is stabilized when the 
star is spun up and the effective gravity is reduced to a certain point (Kulka- 
rni &i Romanova 2009). The first effect, namely, the dependence of the state 
(stable or unstable) on the accretion rate M, may have important observational 
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Figure 2. Left: Accretion in stable regime. The surface is a constant density 
surface, and the lines are sample magnetic field lines. The magnetic axis of 
the dipole [i is inclined relative to the rotational axis fl at 8 = 15°. Right: 
Accretion in the unstable regime at = 5° (from Romanova, Kulkarni & 
Lovelace 2008). 



consequences: a particular star may transition between these two regimes and 
may thus show intermittency of pulsations. The effect of intermittency has been 
observed in a few millisecond pulsars (e.g. Altamirano et al. 2008). This discov- 
ery greatly changes our understanding of accreting magnetized stars and their 
possible observational properties. 




Figure 3. Left: Accretion to a star with a dipole (0 = 45°) and a 
quadrupole (0d = 30°) magnetic field (of comparable amplitudes near the 
star) misaligned at <!> = 90°. A constant density surface and magnetic field 
lines are shown. Right: Same as in the left panel, but showing density con- 
tours in the equatorial plane (from Long, Romanova & Lovelace 2008). 
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4. Accretion to a star with complex magnetic field 

The stellar magnetic field may be more complex than dipole (e.g., Donati et al. 
2007). As a first step we investigated accretion to a star with mixed dipole and 
quadrupole fields: 

, 3(u-f)f — ii 3D. .a 9 3-D . * A . * 

B ( r ) = j + 4^4 ( 5 (° • r) 2 " !) f " ^4 (D • *)D, 

where is the dipole moment, D is the quadrupole moment, and r and D are 
the unit vectors for the position and the quadrupole moment respectively. In 
general, the dipole and quadrupole moments /i and D are misaligned relative 
to the rotational axis fi, at angles O and @d respectively. In addition, they 
can be in different meridional planes with an angle $ between the Q — fi and 
fi — D planes. First, we investigated the case when the moments are aligned 
with each other but not with the rotational axis. We found that in this case, a 
significant amount of matter may flow through the "quadrupole belt" forming 
a ring-shaped spot on the surface of the star (Long et al. 2007). In the more 
general case when the dipole and quadrupole moments are not in the same plane, 
the field is more complicated with a number of poles of different polarity on the 
surface of the star (Long et al. 2008). The accreting matter chooses the most 
energetically favorable path, due to which the funnel streams are often quite 
close to the equatorial plane (see Fig. 3). Work on accretion to a star with a 
higher multipolar field is in progress. 
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